{
 "cells": [
  {
   "cell_type": "code",
   "id": "initial_id",
   "metadata": {
    "collapsed": true,
    "ExecuteTime": {
     "end_time": "2025-04-03T13:58:09.573190Z",
     "start_time": "2025-04-03T13:58:08.442454Z"
    }
   },
   "source": [
    "print(\"begin\")\n",
    "# 用于处理数据\n",
    "import matplotlib\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "from numpy.fft import fft\n",
    "from scipy.signal import butter, lfilter\n",
    "import os,glob\n",
    "\n",
    "# 读取txt数据，每行只有一个数据，并除去里面为0或者格式为2000-0-0 0:5:56.100的行\n",
    "\n",
    "import matplotlib.pyplot as plt\n",
    "import scipy\n",
    "from scipy import signal\n",
    "from scipy.ndimage import median_filter\n",
    "\n",
    "\n",
    "matplotlib.use('TkAgg')\n",
    "x_coords = None\n",
    "# 电刺激开始标识\n",
    "SEMGbeginline = []\n",
    "IMUBeginline = []\n",
    "Tourchbeginline = []\n",
    "\n",
    "SEMG_time = []\n",
    "tourch_time = []\n",
    "window_size = 40\n",
    "step_size = 5\n",
    "\n",
    "def butter_notch_filter(data, fs, fc, bw):\n",
    "    \"\"\"\n",
    "    Apply a notch filter to the data with the given center frequency, bandwidth, and sample rate.\n",
    "\n",
    "    :param data: The input data to be filtered.\n",
    "    :param fs: The sample rate of the input data.\n",
    "    :param fc: The center frequency of the notch filter.\n",
    "    :param bw: The bandwidth of the notch filter.\n",
    "    :return: The filtered data.\n",
    "    \"\"\"\n",
    "    # Calculate the notch filter coefficients\n",
    "    b, a = butter(4, [fc - bw / 2, fc + bw / 2], fs=fs, btype='bandstop')\n",
    "\n",
    "    # Apply the notch filter to the data\n",
    "    filtered_data = lfilter(b, a, data)\n",
    "\n",
    "    return filtered_data\n",
    "def AMPD(data):\n",
    "    \"\"\"\n",
    "    实现AMPD算法\n",
    "    :param data: 1-D numpy.ndarray \n",
    "    :return: 波峰所在索引值的列表\n",
    "    \"\"\"\n",
    "    p_data = np.zeros_like(data, dtype=np.int32)\n",
    "    count = data.shape[0]\n",
    "    arr_rowsum = []\n",
    "    for k in range(1, count // 2 + 1):\n",
    "        row_sum = 0\n",
    "        for i in range(k, count - k):\n",
    "            if data[i] > data[i - k] and data[i] > data[i + k]:\n",
    "                row_sum -= 1\n",
    "        arr_rowsum.append(row_sum)\n",
    "    min_index = np.argmin(arr_rowsum)\n",
    "    max_window_length = min_index\n",
    "    for k in range(1, max_window_length + 1):\n",
    "        for i in range(k, count - k):\n",
    "            if data[i] > data[i - k] and data[i] > data[i + k]:\n",
    "                p_data[i] += 1\n",
    "    return np.where(p_data == max_window_length)[0]\n",
    "\n",
    "\n",
    "# 对信号进行FFT,横坐标为实际频率\n",
    "def fftsignal(signal):\n",
    "    fft_signal = fft(signal)\n",
    "    freqs = np.fft.fftfreq(len(signal), 1 / 1000)\n",
    "    return freqs, fft_signal\n",
    "\n",
    "def read_IMUdata(file_path):\n",
    "    global IMUBeginline\n",
    "    with open(file_path, 'r') as file:\n",
    "        lines = file.readlines()\n",
    "\n",
    "    IMU1_AngleZ_data = []\n",
    "    IMU2_AngleZ_data = []\n",
    "\n",
    "    for line in lines[1:]:\n",
    "        line = line.strip()\n",
    "        if not line:\n",
    "            continue\n",
    "        elif line == \"FESBegin\":\n",
    "            IMUBeginline.append(len(IMU1_AngleZ_data))\n",
    "        else:  # 正常数据识别\n",
    "            parts = line.split()\n",
    "\n",
    "            if parts[0] == \"IMU_1\":\n",
    "                IMU1_AngleZ_data.append(float(parts[5]))\n",
    "            elif parts[0] == \"IMU_Sta\":\n",
    "                IMU2_AngleZ_data.append(float(parts[5]))\n",
    "\n",
    "    if len(IMU1_AngleZ_data) == 0 :\n",
    "        return IMU2_AngleZ_data\n",
    "    elif len(IMU2_AngleZ_data) == 0:\n",
    "        return IMU1_AngleZ_data\n",
    "    else:\n",
    "        IMU1_AngleZ_data = np.array(IMU1_AngleZ_data) - np.array(IMU2_AngleZ_data)\n",
    "        return IMU1_AngleZ_data\n",
    "\n",
    "def read_Tourchdata(file_path):\n",
    "    global Tourchbeginline, tourch_time\n",
    "    time_Tourch_z = []\n",
    "    Tourch_z = []\n",
    "    \n",
    "    with open(file_path, 'r') as file:\n",
    "        lines = file.readlines()\n",
    "    for line in lines[1:]:\n",
    "        line = line.strip()\n",
    "        if not line:\n",
    "            continue\n",
    "        elif line == \"FESBegin\":\n",
    "            IMUBeginline.append(len(Tourch_z))\n",
    "        else:  # 正常数据识别\n",
    "            parts = line.split()\n",
    "            time_Tourch_z.append(float(parts[0]))\n",
    "            Tourch_z.append(float(parts[1]))\n",
    "    Tourch_z = np.array(Tourch_z)\n",
    "    Tourch_z = Tourch_z*9.8 \n",
    "    Tourch_z.round(5)\n",
    "    # Tourch_z.astype(float)\n",
    "    tourch_time = np.array(time_Tourch_z)\n",
    "    return Tourch_z\n",
    "    \n",
    "    \n",
    "    \n",
    "\n",
    "def read_SEMGdata(file_path):\n",
    "    global SEMGbeginline, SEMG_time\n",
    "    data = np.loadtxt('实验记录SEMG20241213190046.txt')\n",
    "    #读取n*4矩阵\n",
    "    return data\n",
    "    \n",
    "        \n",
    "\n",
    "\n",
    "# 对读取到的数据进行加窗，每个窗口有50个数据，并计算RMS，每次移动step_size个数据\n",
    "def windowing_and_rms(data, window_size, step_size):\n",
    "    \"\"\"\n",
    "    对读取到的数据进行加窗，每个窗口有window_size个数据，并计算RMS，每次移动step_size个数据\n",
    "    :param data:\n",
    "    :param window_size:\n",
    "    :param step_size:\n",
    "    :return:\n",
    "    \"\"\"\n",
    "    windowed_data = []\n",
    "    global x_coords\n",
    "\n",
    "    for i in range(0, len(data) - window_size, step_size):\n",
    "        window = np.array(data[i:i + window_size])\n",
    "        window_rms = np.sqrt(np.mean(window ** 2))\n",
    "        # 对window_rms进行中位滤波\n",
    "        # window_rms = median_filter(window_rms, 10)\n",
    "        windowed_data.append(window_rms)\n",
    "\n",
    "    # 放缩data的x坐标w使之与windowed_data的x坐标对应\n",
    "    x_coords = np.arange(len(windowed_data)) * step_size\n",
    "\n",
    "    return windowed_data\n",
    "# 50Hz工频陷波器\n",
    "def filter_50Hz(data):\n",
    "    # 50Hz陷波器\n",
    "    b, a = signal.butter(4, [50, 100] / (250 / 2), btype='bandstop')\n",
    "    filtered_data = signal.filtfilt(b, a, data)\n",
    "    return filtered_data\n",
    "\n",
    "\n",
    "# 提取包络线特征\n",
    "def envelope_feature(data, window=20):\n",
    "    envelope = []\n",
    "    data = np.array(data) - 1620\n",
    "    # 每window个数据计算均值并提取包络特征\n",
    "    for i in range(0, len(data) - window, window):\n",
    "        datawindow = data[i:i + window]\n",
    "        envelope.append(np.abs(np.mean(datawindow)))\n",
    "    return envelope\n",
    "\n",
    "\n",
    "def plot_datasimple2(data, windowing_and_rms):\n",
    "    fig, axs = plt.subplots(2, 1, figsize=(10, 8))\n",
    "    axs[0].plot(data)\n",
    "    axs[0].set_title('Original Data')\n",
    "    # x_coords,\n",
    "    axs[1].plot(windowing_and_rms)\n",
    "    axs[1].set_title('RMS Data,winSize = {},stepSize = {}'.format(window_size, step_size))\n",
    "\n",
    "    plt.tight_layout()\n",
    "    plt.show()\n",
    "\n",
    "\n",
    "\n",
    "# 参数化绘制subplot的函数\n",
    "\n",
    "\n",
    "\n",
    "# freqs, fft_signal = fftsignal(SEMGdata)\n",
    "# ... 其他特征\n",
    "plt.rcParams['font.sans-serif'] = ['SimHei']  # 用来正常显示中文标签\n",
    "plt.rcParams['axes.unicode_minus'] = False  # 用来正常显示负号]\n",
    "\n",
    "\n",
    "semgdata = read_SEMGdata('实验记录SEMG20241213190046.txt')\n",
    "\n",
    "print(semgdata.shape)"
   ],
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "begin\n",
      "(10010, 4)\n"
     ]
    }
   ],
   "execution_count": 1
  },
  {
   "metadata": {
    "ExecuteTime": {
     "end_time": "2025-04-03T13:58:16.119342Z",
     "start_time": "2025-04-03T13:58:11.659051Z"
    }
   },
   "cell_type": "code",
   "source": [
    "def plotrawsEMG(data):\n",
    "    plt.plot(data)\n",
    "    plt.xlabel('Time')\n",
    "    plt.ylabel('EMG Amplitude')\n",
    "    plt.title('Raw EMG Signal')\n",
    "    plt.grid(True)\n",
    "    plt.show()\n",
    "plotrawsEMG(semgdata[:,1])\n",
    "def plot_sEMG_data_4channel(data):\n",
    "    fig, axs = plt.subplots(4, 1, figsize=(10, 8))\n",
    "    for i in range(4):\n",
    "        axs[i].plot(data[:, i])\n",
    "        axs[i].set_title('Channel {}'.format(i+1))\n",
    "        axs[i].grid(True)\n",
    "        axs[i].set_xlabel('Time')\n",
    "        axs[i].set_ylabel('EMG Amplitude')\n",
    "    fig.tight_layout()\n",
    "    plt.show()\n",
    "# plot_sEMG_data_4channel(semgdata)"
   ],
   "id": "932884dd6394faed",
   "outputs": [],
   "execution_count": 2
  },
  {
   "metadata": {
    "ExecuteTime": {
     "end_time": "2024-12-16T11:49:22.547678Z",
     "start_time": "2024-12-16T11:49:02.017715Z"
    }
   },
   "cell_type": "code",
   "source": [
    "from mpl_interactions import interactive_plot\n",
    "\n",
    "meannum = np.mean(semgdata[:,1])\n",
    "print(meannum)\n",
    "singalchsemg = semgdata[:,1] - meannum\n",
    "\n",
    "plotrawsEMG(singalchsemg)\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "def fftsignal(signal):\n",
    "    \"\"\"\n",
    "    对信号进行FFT,横坐标为实际频率\n",
    "\n",
    "    :param signal: 输入信号\n",
    "    :return: 频率数组, FFT信号数组\n",
    "    \"\"\"\n",
    "    fft_signal = fft(signal)\n",
    "    freqs = np.fft.fftfreq(len(signal), 1 / 1000)  # 假设采样率为1000 Hz\n",
    "    print(freqs)\n",
    "    return freqs, fft_signal\n",
    "# 假设 singalchsemg 已经被定义并加载\n",
    "# singalchsemg = semgdata[:, 1]\n",
    "\n",
    "\n",
    "# 对信号进行FFT变换\n",
    "freqs, fft_signal = fftsignal(singalchsemg)\n",
    "\n",
    "\n",
    "\n",
    "def plot_fft_signal(freqs, fft_signal):\n",
    "    \"\"\"\n",
    "    绘制FFT信号\n",
    "\n",
    "    :param freqs: 频率数组\n",
    "    :param fft_signal: FFT信号数组\n",
    "    :return: None\n",
    "    \"\"\"\n",
    "    # 绘制FFT结果\n",
    "    fig, ax = plt.subplots(figsize=(10, 6))\n",
    "    plt.plot(freqs, np.abs(fft_signal))\n",
    "    \n",
    "    plt.title('FFT of SEMG Signal')\n",
    "    plt.xlabel('Frequency (Hz)')\n",
    "    plt.ylabel('Amplitude')\n",
    "    plt.xlim(0,500)  # 通常只显示低频部分\n",
    "    plt.grid(True)\n",
    "    \n",
    "    plt.show()\n",
    "\n",
    "# 20Hz高通滤波器\n",
    "b, a = signal.butter(4, 100 / (1000 / 2), btype='highpass')\n",
    "filtered_signal = signal.filtfilt(b, a, singalchsemg)\n",
    "# 300Hz低通滤波器\n",
    "b, a = signal.butter(4, 400 / (1000 / 2), btype='lowpass')\n",
    "filtered_signal = signal.filtfilt(b, a, filtered_signal)\n",
    "\n",
    "plotrawsEMG(filtered_signal)"
   ],
   "id": "7765e9ae4e3ada58",
   "outputs": [],
   "execution_count": null
  },
  {
   "metadata": {
    "ExecuteTime": {
     "end_time": "2024-12-16T12:09:13.522291Z",
     "start_time": "2024-12-16T11:52:19.421416Z"
    }
   },
   "cell_type": "code",
   "source": [
    "\n",
    "\n",
    "def compare_fft_signal(original_signal, filtered_signal):\n",
    "    \"\"\"\n",
    "    比较原始信号和滤波信号的FFT变换结果\n",
    "\n",
    "    :param original_signal: 原始信号\n",
    "    :param filtered_signal: 滤波信号\n",
    "    :return: None\n",
    "    \"\"\"\n",
    "    # 对原始信号和滤波信号进行FFT变换\n",
    "    fig, axs = plt.subplots(2, 1, figsize=(10, 8))\n",
    "    freqs, fft_original = fftsignal(original_signal)\n",
    "    freqs, fft_filtered = fftsignal(filtered_signal)\n",
    "    axs[0].plot(freqs, np.abs(fft_original))\n",
    "    axs[0].set_title('Original Signal FFT')\n",
    "    axs[1].plot(freqs, np.abs(fft_filtered))\n",
    "    axs[1].set_title('Filtered Signal FFT')\n",
    "    axs[0].set_xlim(0, 500)\n",
    "    axs[1].set_xlim(0, 500)\n",
    "    fig.suptitle('FFT Comparison')\n",
    "    plt.xlabel('Frequency (Hz)')\n",
    "    plt.ylabel('Amplitude')\n",
    "    plt.tight_layout()\n",
    "    plt.show()\n",
    "compare_fft_signal(singalchsemg, filtered_signal)\n",
    "\n",
    "\n",
    "def filter_signal_compare(origin, newsignal):\n",
    "    \"\"\"\n",
    "    对信号进行滤波，并比较原始信号和滤波信号\n",
    "\n",
    "    :param signal: 输入信号\n",
    "    :param lowcut: 低通滤波截止频率\n",
    "    :param highcut: 高通滤波截止频率\n",
    "    :param order: 滤波器阶数，默认为5\n",
    "    :return: None\n",
    "    \"\"\"\n",
    "    # 原始信号\n",
    "    plt.subplot(2, 1, 1)\n",
    "    plt.plot(origin)\n",
    "    plt.title('Original Signal')\n",
    "    plt.xlabel('Time')\n",
    "    plt.ylabel('Amplitude')\n",
    "    plt.grid(True)\n",
    "    # 滤波信号\n",
    "    plt.subplot(2, 1, 2)\n",
    "    plt.plot(newsignal)\n",
    "    plt.title('Filtered Signal')\n",
    "    plt.xlabel('Time')\n",
    "    plt.ylabel('Amplitude')\n",
    "    plt.grid(True)\n",
    "    plt.tight_layout()\n",
    "    plt.show()\n",
    "    plotrawsEMG(newsignal)\n",
    "    plotrawsEMG(origin)\n",
    "filter_signal_compare(singalchsemg, filtered_signal)"
   ],
   "id": "4d90edabebd416dc",
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[ 0.         0.0999001  0.1998002 ... -0.2997003 -0.1998002 -0.0999001]\n",
      "[ 0.         0.0999001  0.1998002 ... -0.2997003 -0.1998002 -0.0999001]\n"
     ]
    }
   ],
   "execution_count": 47
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}
